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Recently, three-dimensional Stirling engine simulations have been ac- 
complished utilizing commercial Computational Fluid Dynamics software. 
The validations reported can be somewhat inconclusive due to the lack 
of precise time accurate experimental results from engines, export con- 
trol/proprietary concerns, and the lack of variation in the methods uti- 
lized. The last issue may be addressed by solving the same flow problem 
with alternate methods. In this work, a comprehensive examination of the 
methods utilized in the commercial codes is compared with more recently 
developed high-order methods. Specifically, Lele’s Compact scheme and 
Dyson’s Ultra Hi-Fi method will be compared with the SIMPLE and PISO 
methods currently employed in CFD-ACE, FLUENT, CFX, and STAR- 
CD (all commercial codes which can in theory solve a three-dimensional 
Stirling model although sliding interfaces and their moving grids limit the 
effective time accuracy). We will initially look at one-dimensional flows 
since the current standard practice is to design and optimize Stirling en- 
gines with empirically corrected friction and heat transfer coefficients in 
an overall one-dimensional model. This comparison provides an idea of the 
range in which commercial CFD software for modeling Stirling engines may 
be expected to provide accurate results. In addition, this work provides a 
framework for improving current one-dimensional analysis codes. 
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Nomenclature 


/ 3 Nondimensional Wavenumber 

e Turbulent Dissipation Rate, (m 2 /s 3 ) 

r) Kolmogorov Length Scale, (m) 

l_i Dynamic (Molecular) Viscosity, (Ns ■ m~ 2 ) 

v Kinematic Viscosity, (// ■ p -1 ) 

p Density, (kg- m ~ 3 ) 

t w Surface Shear Stress, (kg ■ m -1 s -2 ) 

c Constant Convective Velocity 

c4o0 UHF Method - 4 point stencil - no derivatives 

c4ol UHF Method - 4 point stencil - one derivative 

c4o2 UHF Method - 4 point stencil - two derivatives 

c4o3 UHF Method - 4 point stencil - three derivatives 

C p Specific Heat Constant Pressure 

C v Specific Heat Constant Volume 

G Complex Amplification Factor 

k Wavenumber 

L Length of Numerical Domain 

r Courant Number 

u T Friction Velocity, ( kg ■ m~ 3 ) 

v Von Neuman Number 

CFL Courant-Friedrichs-Lewy 


I. Introduction 

Power conversion with free-piston Stirling engines 1 promises to deliver high efficiency, low 
mass solutions for longer and more varied space missions. 2 In addition to using advanced 
high-temperature materials to increase the Carnot temperature ratio, it is anticipated that 
advanced computational fluid dynamics (CFD) will help to identify the following losses 4-6 
(also shown in figure 1): 

1. Inefficient heat exchange and pressure loss in the heat exchangers (heater, regenerator, 
and cooler) 

2. Gas spring and working space loss due to hysterisis and turbulence, 

3. Appendix gap losses due to pumping and shuttle effects, 
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4. Mixing gas losses from unequal temperature distributions or losses from mixing of gas 
streams, or elements of gas at different temperatures, 

5. Conduction losses from the hot to cold regions 

In addition, the following artificial numerical losses must be considered when computa- 
tional simulations are performed (also shown in figure 1): 

1. Moving/deforming mesh losses from repeated low order flow held interpolations, 

2. Transient /Unsteady heat transfer and how loss from inconsistent and inaccurate time 
discretization, 

3. Flow loss from low order approaches resulting in effectively adding artihcial dissipation 
terms along sliding interfaces, at structured/unstructured grid interfaces, and within 
interior. 

Minimizing those artihcial losses is best accomplished through higher order approaches. 7 
While this approach is common in aeroacoustics, computational electromagnetics, and exte- 
rior how problems, high order techniques have not yet been applied to simulating a Stirling 
device. Moreover, the following difficulties are often encountered when using high-order 
approaches: 

1. Generation of high-order, smooth, body-htted grids around complex configurations can 
be difficult. 8 

2. High-order formulations can lack nonlinear robustness. 8 

3. The general usefulness of high-order methods is limited by hrst order accurate shock 
capturing. 9 

Fortunately, with the exception of the possibly random geometry in the regenerator, the 
free-piston design is essentially smooth and admits curvilinear structured grids (with some 
geometry simplification). The issue of nonlinear robustness (i.e. maintaining design accu- 
racy with nonlinear equations) is an open issue, but preliminary results are encouraging. 
And finally, the working gas is subsonic and shockless throughout the entire region 10 (How- 
ever, steep temperature gradients can exist at solid/fluid interfaces). For these reasons, a 
high-order approach is investigated for ’’whole engine” Stirling analysis and compared to 
commonly used techniques. 
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Schematic with Springs and Dampers 


Effect of Heating & Cooling 


Artificial Transient Heat Loss 



/ Bearing System 

Moving Mesh Loss 



Working Fluid (bounce space and working space) 


Figure 1. Schematic of Actual and Artificial Numerical Losses in a Free-Piston Convertor 


II. Description of the Problem 

The dual opposed configuration shown in figure 2 n> 12 is being developed for multimission 
1 

uses. 





Figure 2. Dual Opposed Stirling Convertors Reduce Vibration 

Many methods in general use stop at A th order accuracy for time dependent problems 
since they use Runge-Kutta methods, ffigh-order Runge-Kutta methods become notoriously 
difficult to derive because the number of nonlinear order conditions that need to be solved 
grows exponentially (i.e. , a 12 th order method has 7813 nonlinear order conditions). The 
advantages of using Runge-Kutta methods at orders less than 6 are commonly cited as 
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flexibility, large stability limits, and ease of programming . 14 The practical limit on their 
order has been an impediment to the analysis of their use in high order approaches for time 
dependent applications . 15 

In this paper we use a series of explicit, local, high order methods which have the same 
order of accuracy in space as in time 16, 1 ' for inviscid flow (lower accuracy in time for viscous 
flows). These methods use Hermite interpolation on stencils that are four points wide, and 
a Cauchy-Kowalewski recursive procedure 18 for obtaining time derivatives from the space 
derivatives of the interpolant. The time derivatives are then used to advance the primitive 
variables and their spatial derivatives in time with a Taylor series expansion. This general ap- 
proach is called the Modified Expansion Solution Approximation (MESA) method 19 and the 
new finite volume variation of this is called the Ultra HI-FI (UHF) method . 20 This method 
can be used to derive and implement algorithms with arbitrarily high orders of accuracy 
in multiple space dimensions if their complexity is properly managed and the computer’s 
floating point precision is sufficiently high . 21 

First, some of the known exact solutions of the viscous Burger’s equation are provided and 
the linear case is solved with both state-of-practice Compact schemes, new UHF methods, 
and various commercial code solvers. The one-dimensional Navier-Stokes equations reduce 
to the linear viscous Burger’s equation in certain circumstances and provide a means for 
testing both heat transfer effectiveness and turbulent transition efficiency of each method. 

III. Exact Solutions For Method Comparison 

Since the nonlinear Navier-Stokes equations generally do not have exact solutions nor 
known stability limits, preliminary development and testing of new numerical methods is best 
accomplished by starting with the viscous Burger’s equation. This equation describes flow 
behaviour in specialized circumstances, but more importantly, its mathematical properties 
are very similar to the full Navier-Stokes equations and it admits exact solutions. 

For future reference and convenience, some of the known exact solutions are shown below 
(only the linear viscous Burger’s equation will be required in this work). 

A. Complete Nonlinear Viscous Burger’s Equation 

The viscous Burger’s equation is written as: 

du du d 2 u 
dt U dx ^ dx 2 

where u is the convective velocity term and /i can be considered the dynamic viscosity. 

Exact steady-state solution, lim t ^oou(x,t) exists for the case with boundary conditions: 
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w(0, t) = uo 
u(L } t) = 0 


and it is given by: 


1 — expluRe^x / L — 1)] 

0 1 + exp[uRe l(x / L — 1)] 

where ReL = (note this is a modified Reynold’s number) 
u is a solution of the equation 


( 2 ) 

(3) 

(4) 


— — = exp (~uRe L ) (5) 

u + 1 

Other solutions for the nonlinear viscous burgers equation are: 

1 . 

U t + uu x - nu xx = 0,fjL = 0.1, X e (0, 1) (6) 


u(x, 0) = 0, w(l, t) = 


-tank f — j ,u(0,t) = 0 


has exact solution 


x 


u(x,t) = —tank — 
y z/i 


2. Fully nonlinear equation: 


(7) 


( 8 ) 


u t + uu x - - ( uu x ) x = 0, x G (0, 1) 


with initial and boundary conditions as: 


u(x, 0) = exp x ’, w(0, t) = 1, u(l, t) = e 


(9) 


( 10 ) 


has the following exact solution u(x,t ) = exp 2 ’. 

3. An exact solution of the nondimensional form of the Nonlinear Burger’s equationn 22 


du du d 2 3 u 
— — b u — — = / 1 


dt dx ^ dx 2 


( 11 ) 
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Nondimensionalized by: 


producing: 



uL a 

— ’** = J2 
li L 2 


t 


u* „ du* d 2 u* 

l- u * — 

t* dx* dx* 2 


One stationary solution is: 


2 sinh x* 
cosh x* — exp - ** 


( 12 ) 


(13) 


(14) 


B. Linear Viscous Burger’s Equation 

For simplicity and more thorough stability analyses, the viscous Burger’s equation may be 
linearized with constant convective velocity, c, and dynamic viscosity, [i: 


du du d 2 u 
dt ° dx ^ dx 2 


(15) 


The exact steady-state solution with the same boundary conditions as in Eqs.(2) and (3) 


is: 


U = Uq 


1 - exp [R l (x/L - 1)] 
1 - exp (~R l ) 


where Rl = ^ (modified Reynolds number). 

The exact solution for the linearized equation with initial condition, u(x, 0) 
and periodic boundary conditions is: 


(16) 


sin(kx), 


u(x, t ) = exp {—k 2 /it) sin k{x — ct ) 


(17) 


This is useful for evaluating the temporal accuracy of a method and this will be used in 
comparing Compact and UHF techniques. 23 

In addition, this equation form also describes the time-accurate temperature distribution 
in a moving solid or within a moving fluid in a channel in which case fi — a — is 
interpreted as the thermal diffusivity, and T is the temperature: 


dT dT _ d 2 T 
dt U dx a dx 2 


(18) 


This will be used for comparing the heat transfer capabilities of Compact, UHF, and 
commercial code solvers. 
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C. 2-D Nonlinear Viscous Burgers’ Equation 

The 2-D Burger’s equation with nonlinear convection terms, 24 


with 


du du du 
m +u dx +l % 


d 2 u 
dx 2 



- fx 


0 


dv dv dv 
dt +U dx +V 8y 


,d 2 v 

'dx 2 


6 ' 2 _v 
dy 2 


~fy 


0 


has exact solution: 


fx = 

fv = 


1 


1 

(uf + 


x 2 + 2 xy 

(1 + 1 ) 

y 2 + 2 xy 


(1 + 1) 


+ 3 x 3 y 2 — 2 vy 
+ 3 y 3 x 2 — 2vx 


(19) 

( 20 ) 


( 21 ) 

( 22 ) 


u = 


1 

1 + 1 


V = 


1 

1 + t 


+ x 2 y 

(23) 

+ xy 2 

(24) 


D. 2-D Linear Burger’s Equation 

The 2D (Burger’s) linear convection and diffusion equation: 25 


u t + c{u x + Uy) - n{u xx + Uyy ) = 0, (x , ?/) G (~1, l) x (~1, l); c = 1, n = 0.01 (25) 


with the initial condition u(x, y, 0) = sin(7r(x + y)) and periodic boundary condition has 
the exact solution is: 

u(x, y, t ) = exp _27r/it sin(7r(x + y — 2 ct)) (26) 

Additional exact solutions may be found in the paper by Benton and Platzman. 26 


IV. Application of Compact Scheme 

The currently accepted state-of-the-art approach to high fidelity numerical simulations 
is based on Lele’s 6 th order Compact scheme. 2 ' This approach implicitly solves the spatial 
derivative terms and utilizes standard Runge-Kutta time advance. Each time step requires 
solving a tridiagonal matrix given by the equation: 
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(du\ (du\ (du\ 

“ \dx) + (& J, + “ (aS J i+1 - “ 2Ax + b iKx 

with a = 1/3, a = 14/9, b = 1/9. 

Similarly, the second order derivatives are found by: 


( 27 ) 


a 


'd 2 u 

dx 2 


i— 1 


'd 2 u\ /^uX _ ttj + i - 2uj + ^i +2 - + ttj-2 

dx 2 ) i a \dx 2 ) i+l a Ax 2 4Ax 2 


(28) 


with a = 2/11, a = 12/11, b = 3/11. 

The 4 th order Runge-Kutta method is described by: 24 


R(u) 

CU X + f-^^XX 

(29) 


= u n + ^R n 

(30) 

u< 2 > 

= + 

(31) 


= u n + AtR 2 

(32) 

u n+1 

= u n + — (R n + 2 R {1) + 2 R {2) + i? {3) ) 

6 v / 

(33) 


(34) 


with 

RW = i?(w (1) ) 

R (2) = R(u {2) ) 

R (3) = R(m (3) ) 

The single step 6 th order Compact scheme with A th order Runge-Kutta on a domain 
[—4, 4] is given by: 


u? +1 = 


(35) 


( 


/ 6296179542799261293r 4 , 9881647455875381 Ur 3 _ 31068146188125r 3 , 55178557611u 2 r 2 1 , 

( 2934584609881498112 _l_ 18549839506204160 42158726150464 _l_ 19381094656 J a ~ A ' 

( 287760797433ur 2 , 324338625r 2 , 401096399245ii 3 r \ , 

V 133245025760 2422636832 ' 1_ 583491009024 ) U ~ A 

/ 3588625vV , 573327»r _ 75r _ 1230500161u 4 _ 268897n 3 , 16675u 2 _ 11?; \ , 

( 3898048 4 1093840 8701 24196548096 1714608 4 127008 180 J U “ 4 _l_ 

76938698802699r 4 _ 6211549394313?;?- 3 , 20491256151r 3 _ 186196891t) 2 r 2 , 1236331119 ot 2 _ 1463625r 2 \ 
12527254840352 860382166336 ~ l ” 7918618736 40043584 “ l ” 346090976 2502724 J 


U- 3 + 
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1778984341n 3 r , 77947t> 2 r _ 65875ur , 75 r , 454075<; 4 _ 203259m 3 _ 1025v 2 , 25?; \ , 

2210193216 ~ l ” 44296 69608 1582 3 4148928 24004512 3528 252 J U ~ 3 

/ 12528089829r 4 , 7466191437?;r 3 _ 22282569r 3 , 441v 2 r 2 _ 18974931?;r 2 , 103671r 2 \ , 

( 1124897312 _l_ 646417856 3652264 _l_ 484 2860256 ~ l " 47432 J a ~ 2 _l_ 


111033673^ _ 49^ 2 r , 15479fr _ 75r 

66975552 66 6328 308 


2401d 4 , 267775v 3 , 49?> 2 
39366 889056 ~ l ” 162 



( 


( 


21 1569041 160693r 4 _ 7860884553897w 3 , 85890653865r 3 , 186196891?; 2 ?- 2 , 2047213695?;r 2 _ 13246263r 2 \ 
12527254840352 860382166336 “ l " 7918618736 “ l " 40043584 346090976 2502724 ) U ~ 1 


3480659477ti 3 r _ 77947?; 2 r _ 198763tir , 2637r _ 454075?; 4 _ 9562127?; 3 , 1025i; 2 , 221?; \ , 

2210193216 44296 69608 _l_ 1582 4148928 24004512 ' 1_ 3528 252 ) LL ~ 1 ' 


( 28605737830234834755r 4 _ 278097527037915r 3 _ 72837778155 v 2 r 2 , 8619732963r 2 \ , 

( 1467292304940749056 21079363075232 9690547328 ~ l ” 1211318416 ) U 0 

( I 19447891?< 2 r _ 50889r , 2706289217?; 4 _ 55091» 2 , i A , 

l” 1 ” 5847072 17402 _r 12098274048 63504 ' l I Uq ~r 


21 1569041 160693r 4 , 7860884553897ur 3 , 85890653865r 3 , 186196891?; 2 ?- 2 _ 2047213695'»r 2 _ 13246263 rA 

12527254840352 860382166336 7918618736 40043584 346090976 2502724 J Ul 


( 3480659477VV _ 77947?; 2 r , 198763?;r , 2637r _ 454075t> 4 , 9562127?; 3 , 1025?; 2 _ 22U;\ , 

V 2210193216 44296 ~ l ” 69608 ' 1_ 1582 4148928 _l_ 24004512 ' 1_ 3528 252 ) Ul ' 

( 12528089829r 4 _ 7466191437fr 3 _ 22282569r 3 , 441-t; 2 r 2 , 18974931nr 2 , 103671r 2 \ , 

( 1124897312 646417856 3652264 _l_ 484 _l_ 2860256 ' 1_ 47432 ) ' 

( I 111033673?' 3 r _ 49u 2 r _ 15479fr _ 75 r _ 2401?; 4 _ 267775?; 3 , 49v 2 , 25?; \ , 

V 66975552 66 6328 308 39366 889056 ’ 1 ” 162 ' l_ 84 ) " l_ 

( 76938698802699r 4 , 6211549394313?;?- 3 , 20491256151r 3 _ 186196891?; 2 ?- 2 _ 1236331 119?;r 2 _ 1463625r 2 \ 
V 12527254840352 " l_ 860382166336 ~ l ” 7918618736 40043584 346090976 2502724 ) 


' 1778984341?; 3 r , 77947?< 2 r , 65875?;r , 75 r , 454075?; 4 , 2032591?; 3 _ 1025^ 2 _ , 

v 2210193216 “ l " 44296 ' 69608 ' 1582 _r 4148928 ' 24004512 3528 252 J u 3 ' 

( 6296179542799261293r 4 _ 98816474558753811ur 3 _ 31068146188125r 3 . 

^ 2934584609881498112 18549839506204160 42158726150464 _l_ 

55178557611f 2 r 2 , 287760797433«r 2 , 324338625r 2 \ , 

19381094656 ~ l " 133245025760 ~ l ” 2422636832 j “ l ” 

401096399245t) 3 r _ 3588625i; 2 r _ 573327 vr _ 75 r _ 1230500161t> 4 , 268897i> 3 , 16675t; 2 , llu \ 

583491009024 3898048 1093840 8701 24196548096 ~ l " 1714608 127008 ' 1_ 180 J 


The linearized viscous Burger’s equation (Eq. 15) is solved with this Compact scheme 
with multiple domain sizes to demonstrate the dependence of the implicitly derived spatial 
derivatives (and the stability limit) on the size of the domain. A Fourier stability analysis is 
performed and the stability of the Compact scheme as a function of Courant (r = ^), Von 
Neumann (v = ^|) numbers, and f3 = A xk is derived as: 


&(G) = 


(36) 


28605737830234834755r 4 _ 278097527037915r 3 _ 72837778155 v 2 r 2 , 8619732963r 2 
1467292304940749056 21079363075232 9690547328 ' 1_ 1211318416 


( 


( 


1944789U> 2 r _ 50889r , 2706289217^ 4 _ 55091^ 2 , 

5847072 17402 “ + " 12098274048 63504 “ l ” 


21 1569041 160693r 4 , 85890653865r 3 , 186196891i; 2 r 2 
6263627420176 “ l " 3959309368 “ l " 20021792 


13246263r 2 \ 
1251362 ) 


cos(/?) + 


77947i; 2 r , 2637r 
22148 ~ l ” 791 


454075H 4 , 1025i) 2 \ 

2074464 ~ l ” 1764 ) 


COS (/3) + 


/ 12528089829r 4 
V 562448656 


22282569r 3 , 44m 2 r 2 , 103671r 2 \ 

1826132 _l_ 242 23716 ) 


cos(2 (5) + 


49 v 2 r 75r 

33 154 


240 li; 4 , 49^\ 
19683 ' 81 ) 


cos(2 (3) + 


76938698802699r 4 


+ 


20491256151r 3 

3959309368 


18619689m 2 r 2 


1463625r 2 A 
1251362 ) 


cos(3/3) + 


6263627420176 


20021792 


77947 v 2 r , 75r , 454075?? 4 
22148 _l_ 791 ' 1_ 2074464 


1025?? 2 \ 
1764 ) 


cos(3 (3) + 


( 6296179542799261293?- 4 
( 1467292304940749056 


31068146188125?- 3 , 55178557611?? 2 ?- 2 

21079363075232 9690547328 


+ 


324338625r 2 

1211318416 


) cos(4/3) 


3588625?? 2 f- _ 150r 
1949024 8701 


1230500161?? 4 , 16675?? 2 \ 
12098274048 ' 63504 ) 


cos(4/3) + 1 


The imaginary part of the amplification factor is: 


3(G) = 


(37) 


( 7860884553897 vr 3 
( 430191083168 


2047213695»r 2 _ 3480659477?? 3 ?- , 198763???- , 9562127?? 3 
173045488 1105096608 “ l ” 34804 ‘ 1_ 12002256 


221v \ 
126 ) 


sin (/3) 



7466191437-»r 3 

323208928 


+ 


1897493 m ?- 2 , 111033673ti 3 r 

1430128 “ l " 33487776 


15479?ir 

3164 


267775?' 3 , 25v \ 
444528 ' 42 ) 


sin (2 /3) 


, / 6211549394313??r 3 _ 1236331119w 2 
~ t " l 430191083168 173045488 


1778984341 ?? 3 ?- , 65875 ???- i 2032591?? 3 
1105096608 34804 ~ l " 12002256 


iff) sin (3/3) 



98816474558753811??r 3 

9274919753102080 


+ 


287760797433?? ?- 2 
66622512880 


401096399245?? 3 r 

291745504512 


573327??r , 268897?? 3 , 
546920 857304 


w) sil ?( 4 ?3) 


The amplification factor is simply, G — y'flGj? + ‘CiiO)-, and the range of stability 
(green area) is shown in Fig. 3. Notice how the stable region changes as the domain size, 
maxi (the number of grid points in each direction), changes. Explicit spatial derivative 
operators do not exhibit this behavior since the stencil size remains constant regardless of 
domain size. 


V. Application of UHF Method 

The MESA and Ultra Hi-Fi Methods are actually a procedure for designing ever more 
accurate numerical methods in which additional information is stored at each cell or grid 
point. For this comparison, the solution variable and up to it’s third spatial derivative will 
be stored at each grid. The notation, c4od, represents an UHF method with a 4 point stencil 
and only the solution variable and up to d derivatives on the grid. The basic procedure has 
been previously published for inviscid problems. 20 

The c4o0 UHF method will use a 4 point interpolation stencil to determine spatial deriva- 
tives as shown in Fig 4. A simple 1 st order Taylor series in time is used: 

u. T +1 =u n + u t (At) (38) 

in which the time derivative is found from the governing equation (Eq. 15) to be: 

dii 

— cUx T (39) 

and the solution variables, u, u x , u xx , are interpolated to the center of the four-point 


11 of 33 



(a) MAXI=2 (b) MAXI=3 



(c) MAXI=4 


Figure 3. Lele Compact Linear Viscous Burgers Equation Stability Range 


12 of 33 








Figure 4. UHF Staggered Grid Diagram 


stencil using: 


< = — (~Ui_3 + 9^_i + 9u i+ i - u i+ 3 ) 


u" = 


— Uj_3 + 27ttj_i) — 27w i+ i + w i+ 3 


«L = 


24Ax 

-- 

2 6 2 


-Ui_3 + Ui _l +U i+ 1 -M i+ 3 


2AaA 


(40) 

(41) 

(42) 


This results in a single time equation: 


u. 


n+( 1 / 2 ) _ _ cAt _ 1 

V 2/r 2 24h 

pAt 9cAt 


16.) UU ‘-i + 


2/i- 


8h 


16 


pAt 9cAt 9 

' W + ^)T + TeJ ““-5 + 


9 \ / u At cAt 

+ 777 ) uu i+l + I 7777“ + 


V 2/A 24 h 


1 

16 


uu i+l 


(43) 


And after two time steps, the original ’’unstaggered” grid has been updated with the 
following: 


u 


n+1 


3 


+ 


/ fiAt cAt 1 
V2 /A ~ 2ih ~ 16 

fiAt 9cAt 9 \ / fiAt cAt 1 

' 2/A + 8 h + 16/ Ui ~ 2 V 2/A ~ 24h _ 16 J + 


uAt 9cAt 9 

+ ^r + i 6 ’ + 


2h 2 
/iAt 9cAt 


9 


2 h 2 
jiAt cAt 


8 h + 16 


pAt cAt 1 

2 h 2 24 h 16J \'2h 2+ 24h~ w) ]Ui + 


cAt 


fiAt 9c At 9 \ ///At 

’ 2/A _ 8h + 16/ V 2/A ~~ 24h 
fiAt 9cAt 9 \ 

"2/A + 8 h + 16/ + 


1 

16 


— 777 U i~ 1 + 
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fiAt 
' 2h? 
fiAt 
2h J 


9cAt | 9 
8 h + 16 


+ 2 


/iAt 9cAt 9 
’ 2/t 2 + 8 h + 16 


/iAt t cAt 1 
2 /i 2 + m ~~ i6 


tq+1 + 


9c At 9 \ //iAt cAf 


/ /jA/ cAt 
V2/^ + 24/i 


8/i ' 16/ V 2/i 2 

l^ 2 
16, 


24/i 16 


Mi+2 


^i+3 


(44) 


As was done with the Compact algorithm, a Fourier stability analysis is completed with 
the real part of amplification factor: 

9fc(G) = (1152) -1 (— 292Ck; 2 + 288r(4r — 5) + 9 (348u 2 + 16(1 — 4r)r + 63^ cos(/3) (45) 
-18 (l2n 2 + 16r(4r - 5) + 9) cos(2^) + (9(1 - 8r) 2 + 4v 2 ) cos(3/?) + 738) 


And the imaginary part of the amplification factor: 


3(G) 


1 

48 


v(— 8r + (8r — 


1) cos (/?) + 5)(sin(2/9) 


26 sin(/3)) 


(46) 


The amplihcation factor, |G| = sJ^G) 2 + 3(G) 2 , is plotted as in Fig. 5 as a function 
of Courant and Von Neumann numbers, and (3 = Axk. The left figure (a) unwraps the 
polar plot in right figure (b). As (3 increases the amplihcation factor simply repeats. It has 
been traditional to use a polar plot due to this property. However, as soon as derivative 
information is stored on the grid, the amplihcation factor continues to vary as (3 > 2n. This 
is due to the ability of these techniques to carry ultra-short wave information. Amplihcation 
factors larger than one for any wavenumber imply it is an unstable method. 



Re 



(b) Polar Diagram 


Figure 5. c4o0 - Linear Viscous Burgers Equation 
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Following this same procedure, but adding additional derivatives on the grid to produce 
methods, c4ol, c4o2, and c4o3 results in the algorithms shown in the appendix. 



(a) r=v=.01 


(b) r=v=.05 



(c) r=v=.10 


(d) r=v=.5 


Figure 6. Ampification/Stability Factor Comparison 


The full equations shown in the appendix are required to perform the Fourier stability 
analysis to enable true method comparisons. The exact amplification factor is also shown 
since we know the exact solution. Notice in Fig. 6 that the various UHF methods more 
closely approach the exact ampification factor as the number of solution derivatives on the 
grid increases. Since the UHF methods can resolve ultra-short waves, the wavenumber range 
in the figure could be extended past k = 2ir. 

The stability region for all the UHF methods is shown in Fig. 7. This gives an indication 
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(a) c4ol (b) c4o2 



(c) c4o3 

Figure 7. UHF Linear Viscous Burgers Equation Stability Range 
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of the allowable time step as the dynamic viscosity and convective terms vary. Notice (see 
definition of Von Neuman number given earlier) the role viscosity plays in reducing time 
step size. Fortunately, the viscosity, fi is a small term generally compared to convection, c. 
The Compact scheme has larger allowable time steps as shown in Fig. 3 for a given grid. 
However, more grid is required for the Compact scheme as shown in the next section and 
despite this apparent advantage, a coarser grid actually results in an effectively larger time 
step for the UHF methods. 

VI. Application of the SIMPLE/PISO Methods 

The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) method 28 and Pres- 
sure Implicit with Splitting of Operators (PISO) method 29,30 are the mainstay of commercial 
fluid dynamics solvers. 



(a) Segregated Process 

Figure 8. Overview of Segregated Solution Technique 

An overview of the SIMPLE method is: 

• Start the iterative process by guessing the pressure field. 

• Use those pressure values to determine the velocity from the momentum equations. 

• Determine a pressure correction such that the continuity equation is satisfied. 

• Find corresponding velocity corrections, and use new pressure and velocity. 

• Repeat this until a velocity field is found that does satisfy continuity. 
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The governing equations are linearized to produce a system of linear equations with one 
equation for each cell in the domain. A point implicit (Gauss-Seidel) linear equation solver is 
used in conjunction with an algebraic multigrid method to solve the resultant scalar system 
of equations. 

The time stepping is first order accurate implicit, and the spatial accuracy is second 
order. 

The PISO algorithm moves the repeated calculations required by SIMPLE inside the 
solution stage of the pressure-correction equation to more closely satisfy the continuity and 
momentum equations. The PISO method takes more time per iteration, but often requires 
fewer iterations, particularly for transient problems as will be demonstrated next. 

Finally, the stability limits of these two approaches are large (limited by the need for 
accuracy) compared to the Compact and UHF methods due to implicit time-stepping. 
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VII. Nonlinear Navier Stokes — ID 


We will now test these techniques on an example of the one-dimensional Navier-Stokes 
equations that reduces to a one-dimensional heat transfer problem governed by the linear vis- 
cous Burger’s equation. This is what is solved in Sage/GLIMPS 32 and HFAST 33 to produce 
reasonable one-dimensional Stirling analysis solutions in industry. 

The one-dimensional Navier-Stokes equations can be written as: 23 
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Sage/Glimps 32 utilize a slightly different form: 
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Notice that Stoke’s stress tensor is replaced by source terms. This also simplifies the 
numerics since only a single second derivative, must be calculated compared to three 
in Eq. 47. We will utilize the standard Navier-Stokes form for comparison with commercial 
software. Extension of these ideas to the GLIMPS form is direct. 


A. Convection and Diffusion 

Few exact solutions exist for the Navier-Stokes equations and therefore validating commercial 
codes is typically done with only approximate experimental information. One special case 
that both has an exact solution and yet includes heat transfer physics relevant to oscillating 
Stirling engines will be shown. 
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First, commercial codes operate either in two or three dimensions. However, it is possible 
to cajole the commercial code into solving a one-dimensional problem by solving the full 
Navier-Stokes equations with an inherently one-dimensional problem such as flow through a 
pipe with an initial temperature ’’shock” as shown in Fig. 9. 

600K 300K 

u=l m/s | -»- | -»- | u=l m/s 

-2 m 0 m 2 m 

Figure 9. Heat Transfer Test 

This problem can be solved with either two or three-dimensional solvers and it reduces to 
the one-dimensional Navier-Stokes Eq. 47. Since the density and velocity are constant, both 
the continuity and momentum equations are satisfied. Only the energy equation actually 
needs to be solved: 


dE t d 

dt dx 


((pC„T + p)u- + 0=0 


(53) 


Using, E t = pC v T , and dropping out terms that are zero, we have the following reduced 
form of the energy equation: 


^ dT „ dT , 3 2 T 

pCv^r = pC v u— = k—r 

at ox ox - 


(54) 


Finally, dividing by pC v and for essentially incompressible flows replacing C v with C p , 
we have the following transport (linear viscous Burger’s equation): 


dT dT _ d 2 T 
dt ^ U dx a dx 2 


(55) 


with the thermal diffusivity a = . 

Furthermore, commercial solvers utilize the dimensional form of the Navier-Stokes, but 
by utlizing the following non-dimensionalization (* = nondimensional quantity): 


rj!* _ T * _ x t r 0 * To 

T , X , t ,W U , OL OL 

T 0 L t 0 L L 2 

and using the following nondimensional boundary conditions: 


(56) 


T*(x*,0) = 2 or 1, xe[-2,2], U = [0,1], 
u* = 1.0, Ax* = .1, At* = .01, a* = .03003 


(57) 
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We have an exact solution, by the separation-of- variables, 34 (as N is large): 


N 


T*(x*,t*) = 1 + 0.5- - ^ 


sin 


7 r 


k= 1 


(2k - 1) 


7t(x* ~ U*t *) 


exp[—a*(2k - l) 2 n 2 t* / L 2 ] 
2k — 1 


( 58 ) 


And finally, by choosing the characteristic constants as: 


r 0 = Is, L — 1 m, T 0 = 300/i (59) 

We have the following problem definition for the commercial solvers: 

x G [—2m, 2m], T G [300/i, 600 A'], t G [Os, Is], u = 1 m/s, 

Ax = .lm, At = ,01s, a = .03003m 2 /s, p = 998.2kg/m 2 , 

C p = 4182 J/ (kg • K), k = 125359 J/(s • K), p = 101325 Pa (60) 


This problem is solved and the results are compared as various techniques are applied to 
the problem. In Fig. 10, the exact solution is shown along with four UHF techniques, a 6 th 
order Compact scheme, a segregated spatially and temporally implicit method (SIMPLE) 
used in Fluent, a coupled spatially and temporally implicit method (PISO) used in Fluent, 
and a segregated central difference with a blended (averaged) Euler and Crank-Nicolson 
technique used in CFD-ACE. Fluent only allows 1 st order accuracy in time when moving 
meshes are applied as when the Stirling engine is simulated. However, CFD-ACE allows for 
up to 2 nd order accuracy in time when the meshes are compressed with a spring analogy. 

The best techniques in this test case are the Compact and c4o3 UHF schemes. The 
commercially used solvers are noticeably less able to model time accurate heat transfer. 
However, the coupled solver should be used when commercial codes are applied to oscillating 
Stirling simulations and when available, 2 nd order time accuracy should be used (as when 
using CFD-ACE). 


VIII. Fidelity and Turbulence Transition 

Modeling turbulence transition is a difficult problem due to the large disparity in both 
spatial and temporal scales caused when velocity gradients are high. In the Stirling engine 
velocity gradients are high near walls and regions of sheared flow due to oscillating/reversing 
flows. As the velocity gradients increase, the flow becomes rotational, leading to a vigorous 
stretching of vortex lines, which cannot be supported in two dimensions. 24 For this reason, 
truly turbulent simulations cannot be done in one-dimension. 


21 of 33 


Temperature 0 0 Temperature 0 Temperature 


Exact Solution 

c4o0 

c4o1 

c4o2 

c4o3 

Compact 

Segregated. Imp, Imp 1st 
Coupled. Imp. Imp 1st 
Central. Crank 




Exact Solution 
c4o0 
c4o1 
c4o2 
c4o3 
Compact 
Segregated. Imp. Implst 
Coupled. Imp, Imp 1st 
Central. CranK 

_i i i i i 


(a) Bird’s View 
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(d) Commercial Comparison 
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Figure 10. Heat Diffusion and Convection Comparison 
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The Stirling engine exhibits turbulent and laminar behaviour simultaneously. 35 It is 
desirable therefore to avoid a Reynolds Averaged approach since the time averaged equations 
combined with some turbulence model 36-38 assumes turbulence everywhere. One promising 
approach is Large Eddy Simulation (LES) in which the Navier-Stokes equations are solved 
in time, but with spatial filtering applied, leaving the small eddies still unresolved. Since 
small eddies are essentially isotropic, the modeling is much easier compared to Reynolds 
time averaging. Moreover, the entire flow is bounded by walls making boundary condition 
specification much easier than the typical open domain problems encountered in modeling, 
for example, jet flow turbulence. 

The smallest scales of turbulence are the Kolmogorov scales of length, time and velocity: 39 

n = W 3 M m , T = ( v /eyv (6i) 


where v is the kinematic viscosity, and e is the dissipation rate. The Reynolds number 
in the Kolmogorov region, Re = vrj/u = 1, shows the ratio of inertial and visous forces is 
unity because most of the energy is dissipated in this wavenumber region. 

In small Stirling engines we can estime the smallest scales as follows. With the average 
flow assumptions for Helium in the engine: 

H = 350e — 7(Ns/m 2 ), u = 502e — 6(m 2 /s), k — 278e — 3W/(mK), 
a = 29.9341e — 6m 2 /s, Pr = 0.654, u = 2tt/ = 502.655(rad/s), 

C p = 5.19e3 J/(kgK),T = 700/1, p = 2.6e6 Pa, p = 1.78838(fcg/m 3 ) (62) 

Then from West, 40 for oscillating flow, the average thermal boundary layer thickness, 
\j2a/uj = .345115mm, and the average flow boundary layer thickness, ^2v /u ; = 1.41329mm. 

The maximum surface shear stress may be approximated with this information by (as- 
suming average maximum flow speed of 10 m/s and using the flow boundary layer thickness): 

An 

r w = p — — = ,2476491V/m 2 (63) 

Ay 

The friction velocity, u T = = .372125 m/s, is then used in the following equation for 

dissipation in channel flow: 39 

e « 2w 2 [/ m /. 0253007m) = 54.7325m 2 /s 3 (64) 
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Finally, we can estimate the Kolmogorv spatial wave length as: 


r] = (z^/e) 1 / 4 = 1.23301mm 


( 65 ) 


Assuming a representative engine length of (1 in. = 25.3007mm), the domain consists of 
20 to 30 Kolmogorov wavelengths, or roughly 8000 regions (Kolmogorov boxes) in 3D where 
isotropic turbulence modeling can be employed. This corresponds to a turbulent Reynolds 
number of Re r = u r (.0253007m)/(2z/) = 9.37751. This is low Reynolds number flow and 
appears ideal for applying Large Eddy Simulation. 

One would like to model the smallest turbulent scales, which in the case of LES is the 
Kolmogorov wavenumber range. The allowable time should also be of the same order as the 
Kolmogorov time scale. 


This time scale is not prohibitive, since the time step size is already smaller than this in 
current commercial simulations due to the numerical issues involving moving grids. 

We would like to utilize the most efficient technique to minimize the computational cost. 
We know the minimum wavelengths required to simulate turbulence in the Stirling engine 
and the wave convection/dissipation problem solution from Eq. 17 is used to determine the 
fewest grid points required for each method. The results are shown in Fig. 11. 

This solution represents a traveling (convecting) wave with a dissipating amplitude. Many 
methods will convect at the wrong speed (dispersive error) or will excessively dissipate the 
amplitude (dissipative error). 

A key measure of efficiency is how many grid points are required per wavelength to 
propagate this wave within some predetermined error bound. As k, the wavenumber in 
Eq. 17 , is increased, the number of grid points must increase to properly simulate the wave. 

Currently, Compact schemes are regarded as requiring approximately 6 grid points per 
wavelength for reasonable solutions. A comparison of the Compact scheme and the UHF 
schemes with up to 3 spatial derivatives stored on the grid, as shown in the table (Fig. 11(b), 
demonstrates the Compact scheme is comparable to method c4ol (one solution derivative 
per grid point). This is expected since the Compact scheme also utilizes 1st derivative 
information. However, extending Compact schemes to include higher derivatives involves 
complicated matrices which may be difficult or intractable to solve. 

The c4o3 method can match the results of the Compact scheme using 16 times fewer grid 
points per dimension. Specifically, one additional test of the c4o3 method with k — 8n and 
Ax = .2, the error at t — 1 was 2.44291* - 5, which is more accurate than the Compact scheme 
with considerably fewer grid points per wavelength. Note that the previous comparisons had 


r = 


(z//e) 1/2 = ,00302s 


( 66 ) 
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(a) Fidelity Comparison 


Method 
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.1 
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.4 
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c4ol 
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.4 
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.1 
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c4o2 

.2 
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c4o2 

.4 
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c4o3 

.1 

3.44169 10~ 15 

c4o3 

.2 

2.27818 10" 13 

c4o3 

.4 

3.12925 10 “ n 

Compact 

.1 

1.0993 10~ 6 

Compact 

.2 

7.10929 10" 5 

Compact 

.4 

5.31245 10“ 3 


(b) Errors at t=l 


Figure 11. Wave Propagation Fidelity 
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k = n and here it was multiplied by eight and we use half as many grid points for a net 
change of sixteen grid points per wavelength. 

Clearly, the UHF schemes are more efficient and can be formulated explicitly for easier 
parallelization. This implies smoothly transitioning turbulent flows can be more efficiently 
simulated with UHF techniques. 


IX. Conclusion 

One-dimensional Navier-Stokes equations are currently utilized for Stirling engine design 
and optimization with reasonable success. Recent attempts at multi-dimensional simulations 
have relied upon commercial solvers and this report examined that practice more closely. 

This report has shown that the techniques used in commercial codes for simulating Stir- 
ling engines are not as capable as more recently developed approaches available in the litera- 
ture. Moreover, the unique environment of the Stirling engine in which flow is simultaneously 
turbulent and laminar makes large eddy simulation desirable, while the low Reynold’s num- 
ber, wall bounded flow provides for modest grid requirements and well defined boundary 
conditions. 

Despite the larger stability limit (4 times larger) of Compact schemes for a given grid 
spacing, the UHF method results in an effective time-step that is 4 times larger than Compact 
schemes since the grid can be 16 times coarser per dimension. The 6 th order Compact schemes 
performed well with the heat transfer test and apparently would work well in regions of 
conjugate heat transfer. However, the Compact scheme is not as efficient at predicting 
turbulent transition compared to UHF methods. The c4o3 method performs similarly to 
the Compact scheme for steep temperature gradients (conjugate heat transfer) but is up to 
16 3 = 4096 times more efficient when three-dimensional transitional flows need modeled. 

ft would be desirable to compare c4o4 methods and higher in the future. Future work 
should examine utilizing UHF methods in a steady-harmonic formulation with Detached 
Eddy Simulation and more complicated moving grid tests should be performed. 


26 of 33 


X. Appendix 


In what follows are the full equations used to time advance the solutions in this paper 
and their stability amplification factors. These equations are lengthy, but are provided for 
completeness and to allow for independent verification. 

A. c4ol 

The single step c4ol is given by: 
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The c4ol full step (two staggered half-steps) is: 
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The real part of c4ol amplification factor: 
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(2v — 3) (76n 2 - 9)) sin(2/3)) 

The imaginary part of c4ol amplification factor: 


9(G) = (13824 _1 )(— 81/3 (2560r 3 — 1216(10n + l)r 2 — 8(2n(76n — 99) — 33)r+ (71) 

3(2n + l) (44n 2 -9)) + 

6/3 (9(82n - 39) + 4 (-910n 3 - 3(1104r - 139)n 2 + 5460r(8r - l)n+ 

6r(24r(40r - 23) + 139))) cos (/3) + 

3/3 (23040r 3 + 2112(10n - 3)r 2 - 24 (264n 2 + 38n - 19) r - (2v - 3) (76n 2 - 9)) cos(2/3) 
+2 (l9592n 3 - 972(80r - 3)n 2 - 6(8r(15560r - 2449) + 1845)n+ 
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B 


c4o2 


27(24r(40r(8r - 3) + 9) + 115)) sin(/3) 

+ (2248w 3 + 972(80r - 3)n 2 - 48r(6280r - 281)u - 270n- 
648r(40r(8r - 3) + 9) + 351) sin(2/3)) 


The single step c4o2 (with r = and v = is: 
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The double step from u" to u™ +l can be immediately derived by using two single steps 
above. The equations are ommitted due to length. 
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1. cfo2 Amplification 


The equations for the amplification factor of this method are sufficiently large that their full 
form is not shown. Please see Ref . 41 for the complete form. 
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